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Abstract 

We discuss the computation of the mass of the K and D mesons within the framework of 
A^f = 2 + 1 + 1 twisted mass lattice QCD from a technical point of view. These quantities are 
essential, already at the level of generating gauge configurations, being obvious candidates to 
tune the strange and charm quark masses to their physical values. In particular, we address the 
problems related to the twisted mass flavor and parity symmetry breaking, which arise when 
considering a non-degenerate (c, s) doublet. We propose and verify the consistency of three 
methods to extract the K and D meson masses in this framework. 



1 Introduction 

The framework of maximally twisted mass fermions as an 0{a) improved lattice formulation 
|47| has been proved to be highly successful in recent years. The European Twisted Mass 
Collaboration (ETMC) has adopted this formulation and has carried through a broad research 
program with N{ = 2 flavors of mass-degenerate quarks in various areas of lattice QCD including 
light meson physics [H [21 E] , spectroscopy of light baryons [H [S] , strange and charm physics [UJ 
EE]) -B-physics [HI [10], spectroscopy of static-light mesons [IIl[l2], Isgur-Wise functions [13], 
meson |14 t [T5 l [T6] and nucleon |17j form factors, moments of parton distribution functions |18j . 
neutral [19] and rj' [20] mesons, to — p mass splitting [21], the vacuum polarization tensor [22], 
pion scattering lengths [23], an investigation of the p meson as a resonance [24] or the non- 
perturbative renormalization of quark bilinear operators |25| . 

Particular emphasis has been laid on the cut-off effects appearing at 0{a'^) in the twisted mass 
formulation at maximal twist. These effects have been studied theoretically at tree- level of 
perturbation theory [26], and within the Symanzik approach [321 133|. These analyses suggest 
that isospin breaking effects strongly affect only a limited set of observables, namely the neutral 
pion mass and kinematically related quantities [33]. The same effects have been numerically 
investigated in the quenched approximation [27l [28l [29] , with two dynamical flavors [H [H [5l [SOj 
[31] and with A'^f = 2 + 1 + 1 [38]. All numerical results up to date are in agreement with the 
theoretical conclusions. 

The studies collected so far suggest that the twisted mass formulation at maximal twist is a viable 
realization of QCD on the lattice, with the major advantage of automatic 0{a) improvement of 
physical observables, independently of the specific type of operator considered. Other advantages 
worth to mention are that the twisted mass term acts as an infrared regulator of the theory and 
that mixing patterns in the renormalisation procedure are expected to be simplified. It is 
hence natural to go one step further and include dynamical strange and charm quarks in the 
simulations. The theoretical ground for this has been provided in ref. |34j and first feasibility 
studies have been performed in ref. [35] • In the last years, we have initiated a comprehensive 
research program with dynamical A'f = 2 + 1 + 1 flavors of quarks. Encouraging preliminary 
results were reported in |36[I37]. while a companion paper j38j presents a more detailed analysis 
of the light meson sector for the ensembles used in this paper. 

A difficulty arises in A'f = 2 + 1 + 1 maximally twisted mass lattice QCD when adding a strange 
and a charm quark, due to the explicit violation of the strange and charm flavor quantum 
number conservation. At any non-vanishing value of the lattice spacing, the latter leads to 
the contamination of correlators by unphysical contributions from intermediate states carrying 
the wrong quantum numbers. Moreover, transitions that are not allowed in continuum QCD 
become possible, the consequence being that stable states in the continuum with respect to 
strong interactions, such as the D meson, become resonances. 

In this paper, we provide algorithmic and methodological tools to tackle the problem. In par- 
ticular, we present three techniques, a generalized eigenvalue problem, multiple exponential fits, 
and enforcing parity and flavor symmetry restoration, to compute the physical K and D me- 
son masses. As we will demonstrate below, we find that with all three methods these masses 
can be extracted and results agree among the three methods. The paper is conceived as a 
technical report on these methods, which can in general be applied whenever fiavor symmetry 
breaking occurs. Efforts to implement these techniques in combination with a fiavor diagonal 



Osterwalder-Seiler valence quark action, see e.g. [H [Ml Hll HH] , are ongoing. 

The paper is organized as follows. In section [2] we define the setup, the operators used, and the 
optimization of the correlation matrices. Section [3] describes the determination of the K and D 
meson masses with the three methods. We conclude in section |H 

2 Simulation setup 

2.1 iVf = 2 + 1 + 1 twisted mass lattice QCD 

This work is based on sets of configurations generated by the ETM collaboration |36[ [57] with 
the Iwasaki gauge action |44j and iVf = 2 + 1 + 1 flavors of twisted mass quarks. The light 
degenerate {u, d) quark doublet is described by the standard twisted mass action 05] 



S^Mu[>6'\x^'\U] = a'^^x(^)(x)(Z)w(mo)+i/i75r3)x('H^), (1) 

X 

while for the (c, s) doublet the twisted mass formulation for non-degenerate quarks of |46] has 
been used: 

5F,hcavy[X^"\x^"\C/] = a^ J^ X^'^(^) (^w(mo) + ^^.xTsn + Ts^) X^'^ (^). (2) 

In both cases Dw denotes the standard Wilson Dirac operator 

^w(mo) = ^(7M(v^ + V;)-aV;^V^)+mo, (3) 

while X = (x ;X ) ^^d X = (x i X ) ^'^^ the quark fields in the so-called twisted basis. 
For reasons explained in |35] the same value of the standard quark mass parameter iuq has been 
used in both sectors. 

When tuning the theory to maximal twist, automatic 0{a) improvement for physical quantities 
applies |46[ 147] . This tuning has been done by adjusting mg such that the PCAC quark mass 
in the light quark sector vanishes |38] . 

^mlfrf''' = S T = 0' (4) 



X' 



with the bilinears defined as 

4'^+ = x("^7m75x(') , ^('^+ = X("^75X("^ , P^'^- = X^'^TSX^"). (5) 

At maximal twist, in a massless quark renormalization scheme, the renormalized quark masses 
are related to the bare parameters /ig- and fis by [46] 

m 



« - V(/^.-§w) , m^ = VU + I^ 



Zp^i fJ-a- ^fJ-s] , rn^ = Zp^i i^L^ + ^ ,ns] (6) 



where Zp and Zg are the renormahzation constants of the non-singlet pseudoscalar and scalar 
densities in a massless quark scheme, namely for A'^f = 4 massless Wilson lattice QCD. 

The values of /io- and /X5 have been adjusted in our simulations by requiring that the simulated 
kaon and D meson mass approximately assume their physical values |38j . For this study we 
consider two ensembles, one from each of the currently simulated /3 values, /3 = 1.90 and /3 = 
1.95 [MIEHES], with a light pseudoscalar mass mps ~ 320 MeV in both cases, see Table [TJ 



Ensemble 


/? 


(Ljaf X Tja 


afi 


K 


afia 


af-is 


a 
in fm 


mps 
in MeV 


#of 
gauges 


A40.32 
B35.32 


1.90 
1.95 


32^ X 64 
32^ X 64 


0.0040 
0.0035 


0.163270 
0.161240 


0.150 
0.135 


0.190 
0.170 


0.086 
0.078 


324 
318 


1003 
1042 



Table 1: Summary of the ensembles considered in this paper, more details in \36\ \37\ [38] . 



2.2 Meson creation operators and trial states 

2.2.1 Quantum numbers, physical basis and twisted basis 

We are concerned with computing the mass of the K meson, rriK, and of the D meson, mpi, 
within the setup defined by eqs. ([I|) to (l3|). Both mesons have total angular momentum J = 
and parity V = —. Their quark content is e.g. K^ = ds and D^ = dc. 

Neither heavy flavor nor parity are exact symmetries in A^f = 2 + 1 + 1 twisted mass lattice QCD 
at finite lattice spacing. In particular, the ri-coupling term in eq. ([2]) violates the conservation of 
the strange and charm flavor quantum numbers. Consequently, instead of four different heavy- 
light meson sectors (s, — ), (s,+), (c, — ) and (c, +) there is only a single mixed flavor-parity 
sector (s/c, — /+). Problems arise in particular when one tries to determine mo. In continuum 
QCD the D meson is the lowest state in the (c, — ) sector, while in our setup it is a highly excited 
state in the combined sector {s/c, — /+)■ Notice that, besides the K meson, there are a radially 
excited K state (Ar(1460)), possibly strange mesons with positive parity (^0(800)1 ^oi^^^^)) 
and a number of multi particle states K/Kq + n x vr [¥8]. Hence, for a clean extraction of rriD 
one has to consider sufficiently large correlation matrices, which are able to resolve all these low 
lying states. This is possible in principle. In practice, the separation of the excited states would 
require the determination of correlation matrices with extremely high statistical precision. At 
our currently available statistics, this route seems not to be viable. 

Our approach is instead based on the observation that parity and heavy flavor symmetries 
are restored in the continuum limit, where the twisted mass theory is expected to reproduce 
QCD with N^ = 2 + 1 + 1 quark flavors. In this limit, operators with definite parity |47j 
and flavor quantum numbers projecting onto the physical meson states can be reconstructed (cf. 
section [3^ . As it is shown in the following, these operators can be defined as linear combinations 
of bilinears of the lattice quark fields in the twisted basis. 

In the continuum, or in any chirality preserving lattice formulation [45], the twist transformation 



relating the physical quark fields ■0 and the twisted quark fields x reads 



^{h) 



gJW(75-r3/2 (0 
gii^h75ri/2 (h) 



T^(0g*^i75T3/2 



^W 



y.{h)^iuihl5Tx/2 



(7) 
(8) 



where uji^h are the twist angles in the light and heavy quark sector, respectively. Analogous 
relations hold for operators projecting, in the continuum limit, on trial states with definite 
heavy flavor and parity quantum numbers. In the physical basis, such operators can be chosen 
according tqj 



Oph 



/ 0(3,75) N^ 
ph 

•-^ph 



-zV5W750W 
+-(^('^)^(«) 



(9) 



The twist rotations in eqs. ([7]) and 



/0(S,75)\ 



o 



O 



(c,75) 



o 



{s,l) 



V of '^) 



dS]) relate the twisted basis operators 



to the physical operators of eq. ([9]) as follows 

with the orthogonal twist rotation matrix given by 
/ 



(10) 



(11) 
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cos ^ cos ^ 



sin ^ sin ^ 



sin ^ cos 2 



cos ^ sin 2 



2 
2 
2 
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cos ^ COS ^ 

COS ^ sin ^ 
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cos ^ sin ^ 
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■ sin ^ cos ^ 

- sin ^ sin ^ 

cos ^ cos ^ 



(12) 



/ 



However, when using the Wilson lattice formulation, the operators in eq. (jlOp . with and without 
a 75 matrix, renormalize differently due to the explicit breaking of chiral symmetry. This 
implies that, to be able to build a representation of the chiral group, renormalization factors 
must explicitly be taken into account, and eq. (Ilip only holds for the renormalized counterparts 



'-^ph 



M{ui,L0h)O^ , (O^h)t = {O^)^M^{L0i,ujh), 



(13) 



where the bilinears in eq. (jlOp have been replaced by their renormalized versions, 



O^ = diag(^Zp,Zp,Zs,Zs)0^ 



fZpOi^'^'^\ 
ZpO'^'""^ 
ZsOi''^^ 

\ZsO^^''\ 



(14) 



^For definiteness we identify the light flavor with d. 



and Zp and Zg are the same renormalization factors as in 
7r/2, one has 



). At maximal twist, i.e. w; = w^ 



'-^ph 



R 



( 1 



V 1 



4\ /zpoi^•^^^ 



-1 
1/ 



\Zs0^i^'-\ 



(15) 



A third definition of the quark fields will be useful in the following (where maximal twist applies) , 
obtained by rotating the lattice ^-fields via eqs. (UD and ([SD, where now ui = coh = tt/2. The 
rotated fields would reproduce the physical ones in a theory with exact chiral symmetry and 
Zp = Zs- In the present formulation with broken chiral symmetry, they define instead a "pseudo 
physical basis" (ppb). We denote the rotated fields with ^p'b and introduce the operator 
bilinears in this basis 



Cppb 



'-^ppb 

ppb 

\ '-^ppb 



"*" T^ppb T'^ppb 
_ 7(d) , (c) 
T'^ppb T'^ppb 



/ 1 



V 1 



-1 \ /o?'"^^ 



-1 
1/ 



o 



(c,75) 



o 



(sA) 



\ ot'^ 



(16) 



otherwise written as 



Oppb = A^(^/2,7r/2)0, 



A^mt Ox 



(17) 



The physical operators defined in eq. (J13p . and eq. (jl5|) at maximal twist, project onto states 
that converge to states with definite flavor and parity quantum numbers in the continuum limit. 
Since we aim to determine the ground states of the physical system, in at least two of the four 
sectors, it is appropriate to first build the correlation matrices in terms of the building blocks 
given in eq. (jlOp . We also project to zero momentum by summing over all lattice sites at fixed 
Euclidean time t, 

0?'^Hi) = 77r j;x('')(x,i)rx('^Hx,t), /iG{5,c}, rG{75,l} (77i = ±l,r?^5=±0- (18) 



The corresponding trial states 



\^V\t)) 



o^^^W 1^) 



t, 



enter the correlation matrices 



(19) 



q/..,r.),(^ir,)(*2 - h) = (4'^-^^)(t2)|<^i'^-^^)(ti)) = m{o^^-^^-\t,)) (of 'ri)(to) V) , (20) 

and we introduce the shorthand matrix notation for later use 

C{t2-h) = l0{t2)®{0{h))^). (21) 



Notice also that, due to the discrete symmetries of the twisted mass action in eqs. ([T]) and 
([2]), the correlation C{t2 — ii) is a real and symmetric matrix. Eqs. (jlSp to (j2ip can also 
be generalized to the case of more operators, as for example operators with different levels of 
smearing (see the next section) or Dirac structure. In this case C(i2 — ^i) will be a -D x D matrix 
(D = 4 X n) defined by the larger operator set. An application of this kind will be considered 
in section 13.21 One can easily obtain another set of independent meson creation operators with 
identical quantum numbers by replacing T — > 7or. We found, however, that the corresponding 
trial states have worse overlaps to the low lying states of interest. Therefore, we do not consider 
these operators in the following. To improve the signal-to-noise ratio, we have computed the 
correlators in eq. (j20p by using the one-end trick [U |2] . 

2.2.2 Operator optimization by means of smearing 

To optimize the overlap of the trial states in eq. (|19p with the physical K and D mesons, 
we resort to standard smearing techniques. We use Gaussian smeared quark fields, with APE 
smeared spatial links. Additional details can be found in [11], where the same setup has been 
used. 



We have optimized the smearing by computing effective masses at t = 1 and to = 1 (cf. (I25p ). 
where excited states are suppressed the least, for different values of A'^causs, and KQauss = 0.5, 
Aape = 10, OAPE = 0.5 kept fixed. This optimization is essentially independent on the lattice 
volume and on the light quark mass. Results for /3 = 3.90, L^ x T = 24^ x 48 and /i = 0.0040 
are reported in Figure [TJ Although the suppression of excited states only weakly depends on 
A^Gauss and, therefore, on the width of the corresponding trial states, it is obvious that the D 
meson has a somewhat smaller width than the K meson. Since the D meson is heavier and 
hence more difficult to compute, we focus on optimizing the overlap with the D meson state 
and choose Acauss = 30. An estimate of the corresponding trial state radius R can be obtained 
via [II] 



R I Acauss^Gauss 



a V 1 + 6KGa 



(22) 



yielding Rk ~ 7a ~ 0.60 fm and Rd ~ 5a ~ 0.43 fm (cf. Figure [TJ)). A similar optimization 
for the parameter Aape shows essentially no dependence on the ground state overlap. This is 
exemphfied in Figure [2] corresponding to (3 = 3.90, L^ x T = 24^ x 48 and fi = 0.0100. 

We end up with the following optimized set of smearing parameters for ensemble A40.32: 
Acauss = 30 , KGauss = 0.5 , Aape = 10 , OAPE = 0.5. (23) 

Given the rather mild dependence of the ground state overlap on Acauss and Aape; we use the 
set of parameters in eq. (j23p also for the ensemble B35.32, with only slightly different lattice 
spacing. Sometimes in the following of this paper, we will also consider correlation matrices made 
of local operators, or mixed local and smeared operators. However, the final determination of all 
masses will exclusively be obtained with the correlation matrix made of the smeared operators, 
with the optimized smearing parameters of eq. (j23p . 
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Figure 1: a) The effective masses "T-effective (* ~ -'^'^o = 1) (cf- eq. (p5|) ) for the trial states 
defined in eq. ([H]) as functions of iVcauss for /3 = 3.90, L^ x T = 24^ x 48 and /.t = 0.0040 with 
^Gauss = 0.5, A^APE = 10, aAPE = 0.5. b) The same effective masses as a function of the radius 
of the trial states in lattice units R/a with KQauss = 0.5. 



3 Computation of tjik and ijid 

In contrast to parity and flavor conserving lattice formulations, as the standard Wilson lattice 
QCD, it is not possible to compute correlation functions restricted to a single parity and heavy 
flavor sector in our A'^f = 2 + 1 + 1 twisted mass framework, as outlined in section 12.2.11 While 
the determination of rrix is anyway straightforward, since the kaon is the lowest state in the 
combined heavy flavor and parity sector, the extraction of rriD remains rather problematic, being 
the D meson a highly excited state. Besides computing mx with high precision, we attempt in 
the following to estimate rriD without computing the full low-lying spectrum. We present and 
compare three different methods, all based on the fact that both heavy flavor symmetry and 
parity are only weakly broken, by terms of 0{a). The three methods yield a consistent picture. 



3.1 Method 1: solving a generalized eigenvalue problem 

We consider 4x4 correlation matrices, as defined in eq. (pOl) . computed with the twisted basis 
operators of eq. (jlOp and the optimized smearing parameters given in eq. (|23p . We then solve 
the generalized eigenvalue problem 



J2cjk{ty^\t,to) = J2^Mto)n\tM^^"Ht,to) , t EE t2-ti 



{n), 



(24) 



where k runs over the set (/i, T), h = c,s,T = ±, and obtain the four effective masses "m^^ctn^i 
with n = 0, . . . , 3, by solving 



A(")(^,^o) 
AW(t + l,io) 



e-'"cSLtivc(*.*o)t ^ g-m(j)^^.^^(t,io){T-t) 



,(") 
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e-™effective(*.*0)(t+l) + g-'^cffoctive (*.*o)(T- (t+1)) 



(25) 
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Figure 2: The effective masses ^^effective (* = 1)^0 = 1) as in Figure [U as a function of A'ape for 
/3 = 3.90, L^ X T = 243 X 48 and /jl = 0.0100 with qape = 0.5, iVcauss = 30, KQauss = 0.5. 



with T the temporal extension of the periodic lattice. 

To interpret these effective masses, we assume that heavy flavor and parity breaking effects are 
small. Indeed they are only 0{a), since they originate from the flavor non-diagonal and parity 
odd Wilson term, which is proportional to the lattice spacing. Consequently, for vanishing lattice 
spacing, where heavy flavor and parity are exact symmetries, these correlation matrices would 
be diagonal in the physical basis, because the operators in eq. ([9]) would excite orthogonal trial 
states. Thus, solving the generalized eigenvalue problem as stated in eq. (j24p would directly 
provide the four effective masses with definite heavy flavor and parity. In particular, one of 
them would have associated quantum numbers (c, — ) and would approach a plateau for large 
temporal separation to be identified with the D meson mass. 

At finite lattice spacing in the presence of heavy flavor and parity breaking the four effective 
masses will approach the masses of the four lowest states in the mixed sector {s/c, — /+) for 
large temporal separations. The D meson is not among those states: K and Kq, the radial 
excitations and K/Kq + n x tt states are lighter than the D. At intermediate times, however, 
one of the four effective masses should still be dominated by the D meson and the corresponding 
plateau will give a measure of m£). 

To identify the heavy flavor and parity content of the four effective masses, we flrst note that 
the trial state corresponding to the n-th effective mass is 



l4"^W) = E4"^(t,to)(of(t))V) 



(26) 



When the relations uJi = coh = 7r/2 and Zp/Zs = 1 are approximately fulfllled, one can rotate 
to the pseudo physical basis. By inserting eq. (J16p into the trial state in ()26p and using the 
orthogonality of the twist rotation matrix Mmt at maximal twist of eq. ()17p . yields 



i4"^w) = E G^-t^^"^(*'*o) J^ppb(o tf^) 



t 



(27) 



By sorting the terms in eq. (j27p according to the pseudo physical basis states (C'lpb) 1^)) the 
approximate heavy flavor and parity contents of the trial state corresponding to the n-th effective 



mass can be read off, and it is given by c)^-p^ oc KAImt 'y^"'H*>^o))(/i,r)P- Explicitly, 
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where A^ is a suitable normalization such that 






(32) 



To give a specific example, if c^".,_^ ~ 1, while c^"'.,_^ ~ c}/^^ ~ c,' 



0, the n-th state would 



be interpreted as the D meson. In the continuum limit, where parity and heavy flavor symmetry 

(hT) ~ ^' ^^^ ^^^ others vanishing. 



are restored, each state will have one associated coefficient c, 
Figure [3] shows the first four effective masses "j-eggctive (^ 



0, . . . , 3) as functions of t for the 
ensembles A40.32 (left) and B35.32 (right), while Figure U] shows the approximate heavy flavor 
and parity contents of those states for the ensemble A40.32, measured by the coefficients in 
eqs. (p8|) to (J3T]) . As expected, each one of the effective masses is strongly dominated by and, 
therefore, should correspond to one of the sectors (s,— ), (s,+), (c, — ) and (c, +), which are 
approximately projected by the pseudo physical basis operators associated to the labels (5,75), 
(s,l), (0,75) and (c, 1), respectively. 

To extract the numerical values for uik and rriD, we perform x^ minimizing fits to the corre- 
sponding effective mass plateaus. The fitting intervals [tmin, ^max] sue chosen as follows: 



• ^max = T/2 — 1 = 31 for the K meson. 

• For all the other states tmax is the largest t before which the corresponding effective mass 
is lost in statistical noise (cf. Table [2]). 

• imin is the smallest t fulfilling the following two requirements: 



Z-O ~r J^ — '-min _; '-max- 

All fitting intervals [iminjimax]' with imin 
(x^/dof)max, and we require (x^/dof)max = 



1 < imax < *max, yield a x^/dof < 



2.0. 



By choosing tmm in this way we prevent that effective masses at large t with large statistical 
errors effectively increase the number of degrees of freedom, while not contributing to the 
X^; in practice, the inclusion of these points would allow to fit ranges with too small values 
of imim outside the plateau region. 
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Figure 3: The four effective masses "^effective ^^ functions of t (to = 1) for the ensemble A40.32 
(left) and B35.32 (right). The zoomed in effective masses for the K meson are also shown in 
the bottom graphs. 



Within this method, a systematic error is associated to the determination of the D meson mass, 
due to the fact that the effective mass plateau of the (c, — ) dominated state will finally decay to 
lighter strange states at large times, as a consequence of the heavy flavor and parity breaking. 



We account for this error by taking the difference with a fit in the range \t^ 



l,tr 



and 



we combine statistical and systematic uncertainties in quadrature, where the statistical error is 
obtained by a standard Jackknife analysis. 

The results for mx-, 'mo and the (s, +) state, which for brevity we denote from now on as Kq, 
are collected in Table [2j 

As can also be inferred from Figure [3l we obtain excellent results for rriK. For both ensembles 
the effective mass plateaus extend over more than twenty points, their statistical errors are 
essentially independent of t and the relative errors on mx are ~ 10~^. For rriD the situation 
is more problematic. As shown in Figure [3l the corresponding effective masses are soon lost in 
statistical noise, before they reach unambiguously identifiable plateaus. As mentioned above, 
we add for this a systematic uncertainty. The dominantly (s, +) state does not exhibit a true 
plateau either. One rather observes two different plateaus, and we thus list two results for mx* 
in Table [21 corresponding to two different fitting ranges. A possible explanation might be that 
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Figure 4: Approximate flavor and parity content of the four extracted states as a function of t 
(to = 1) for the ensemble A40.32. Top left: n = 0, mainly (s, — ), i.e. the K meson. Top right: 
n = 1, mainly (s, +). Bottom left: n = 2, mainly (c, — ), i.e. the D meson. Bottom right: n = 3, 
mainly (c, +). The time ranges are the same as for the corresponding effective masses shown in 
Figure [3l 



at small temporal separations t < 10 a positive parity strange meson is seen, while at larger t the 
lighter K + 11 state, with the same strong quantum numbers, dominates. This is also supported 
by the fact that at larger values of the light quark mass a single plateau of rather good quality 
is recovered, see also the results in section] 



3.2 Method 2: fitting the correlation matrix by exponentials 



A complementary approach to determine the heavy-light meson masses is to fit the elements 
of the correlation matrix of eq. (j20p by decomposing them in terms of the eigenstates of the 
Hamiltonian (i.e. the transfer matrix). We consider here the general case with different smearing 
levels, where C{t2 — ti), defined in eq. (j20p . is a D x Z) matrix. When denoting the energy 
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arriK 


t range 


xVdof 


aniK^ 


t range 


xVdof 


arno 


t range 


xVdof 


Ensemble A40.32 


0.2567(2) 


11-31 


0.69 


0.368(32) 
0.473(15) 


14-28 
7-15 


0.92 
1.65 


0.922(11) 


7-14 


0.79 


Ensemble B35.32 


0.2184(3) 


10-31 


0.83 


0.379(28) 
0.446(7) 


12-27 
7-13 


0.54 
1.55 


0.829(8) 


8-16 


0.40 



Table 2: The masses of the K, Kq and D mesons in lattice units obtained by solving a gener- 
alized eigenvalue problem (errors comprise statistical and systematic errors, which are added in 
quadrature). The range and the quality of the fit is also shown. 



eigenstates by |n), n = 1, 2, . . . , M, the matrix elements of C{t2 — ti) can be written as 

M 



Cij{t2-ti) = ^ {i\n)t2 {j\n)t^ 



(33) 



n=l 



with 



(i|n), ^ {n\O^Ht)\n) = (n|(0»(t))V), 



(34) 



where i = 1, . . . ,D labels the operators inserted in the correlation matrix and n = 1, . . . ,M 
counts the eigenstates. Since we consider bosonic operators, we have a periodic time dependence 
on the time extension of the lattice T that can be written as follows 



{i\n)t2 {j\n)ti = (i|n)(j|n)( exp( 



(-(t2 - ti)En) + exp(-(r - t2 + h)En 



(35) 



Here, E^ is the energy of the eigenstate \n) and {i\n) = (i|n)o. In general, the number of 
energy eigenstates is as large as the dimension of the Hilbert space of states. However, for large 
temporal separations t2 — ti,{T — t2 + ti) ^ 1 a few lowest energy states will dominate to a good 
approximation. In this limit, and in analogy with the case of fitting a single correlation function 
with the contributions from a few states, one can fit the matrix of correlation functions with 
the contributions from the set of dominant lowest energy states. In fact, the relevant number of 
energy eigenstates M is small. The number Np of parameters in the fit and the number Nq of 
independent entries of C{t2 — ti) to be fitted are given by 



Np 



M{D + 1) 



Nc 



{tr 



'-min ~r ij 



D{D + l) 



(36) 



where also here imin and imax define the fitting time interval, with (t2 — ti) € [tmin,imax]- The 
minimal set of operators for determining the heavy-light meson masses is given in this case by 
the 4x4 correlation matrix in terms of the operators in eq. (jlOp . The minimal set of states we are 
interested in consists of the K and D mesons. At finite lattice spacing, due to the heavy flavor 
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and parity breaking, the D meson is not stable and does not correspond to an energy eigenstate 
of the lattice theory. However, using the same arguments of section [STTl the D should dominate 
the spectral decomposition of eq. (|33p at intermediate temporal separations. In case of fitting 
the correlation matrix by several exponentials, essential contributions of the lower sectors in the 
D meson channel can be monitored by considering the scalar products of the linear combination 
of the operators obtained from the fit with the rows of the maximal twist matrix. 

Using the pseudo physical basis operators of eq. p^ . one obtains for the coefficients (i|n): 

(i|n) = {n\O^\0)\n) = j; (Mnt)^J^IO^^b(O)ln) . (37) 

j 

Again, assuming that coi = coh ^ '7r/2 and ZpjZg = 1 are approximately verified, the operators 
Oppb should reproduce the physical operators associated to the four channels (s,— ), (s,+), 
(c, — ), (c, +) to a good approximation. In particular, the operator with the same quantum 
numbers of the state |n) should dominate the sum in (j37p . We therefore conclude 

(i|n) ~ Gn(Mrnt) (38) 

to a good approximation, where the proportionality constant G„ is the matrix element of the 
physical operator: 

Gn = mo'^mn). (39) 

It turns out that it is enough to require that the relative signs of the vector components {i\n) 
agree with the signs in the rows of maximal twist rotation matrix A^mt- (A more stringent 
condition on the alignment with the rows of the maximal twist matrix could be imposed by 
requiring the scalar products of the linear combinations of the operators obtained from the fit 
with the rows of the maximal twist matrix to be close to 1, but such a requirement does not 
essentially change the results for the D meson mass.) 

Based on the experience with varying the number of states, we determine the K meson mass 
with a single intermediate state, while good fits for the D meson mass can be obtained for time 
separations around t2 — ti ~ 10 — 12, by using three intermediate states. Taking four states gives 
compatible results, but the signal is lost at smaller distances with consequently larger errors. 
Larger correlation matrices have also been investigated, for instance, 8x8 matrices spanned by 
four Gaussian smeared operators of type (llOp and the corresponding four local operators. In 
this case stable fits with one, three or four states can also be obtained. 

We minimize the uncorrelated x^ 
2 _ v^ f fi{pi,P2,---,PNp)-Xi Y 

where the index i runs over the independent matrix elements to be fitted, Xi and 6Xi are the 
mean value and the error of the matrix element i respectively, and fi{pi,P2, ■ ■ ■ ,PNp) is the 
fitting function depending on Np parameters defined by eqs. (|33p to (j35p . We determined the 
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Figure 5: Masses for the K, (bottom), Kq (middle) and D (top) channels obtained from a 3 x 4 
fit, i.e. a 4 X 4 matrix with 3 states, for the ensemble 535.32. The fit interval is [^1,^2]) with 
values shown by the abscissa. Errors on the masses are also plotted but in most cases are within 
the symbol size - as show by the figure. 



errors of the matrix elements 6Xi and of the fit parameters 6pi by the method in ref. 
Figure [5] illustrates how the extracted masses depend on the fit intervals. We also studied the 
correlated x^ following refs. |5H \52\ 



Nc 



(41) 



where Mj,- = NC-- , with N input data and the estimated covariance matrix 



Ci 



N 



hE(^... 



^i I ( ^i,n — Xj I . 



(42) 



n=l 



It turned out, however, that on our data samples the covariance matrix has a large number of 
almost degenerate tiny eigenvalues of the order of magnitude 10~^®, which cannot be properly 
determined within the present statistical accuracy. The small eigenvalues can be smoothed 
[511I52J . at the price of introducing an uncertainty in the value of Xc- For this reason, we decided 
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to minimize the uncorrelated x ; 3,nd to use the correlated one Xc to estimate systematic errors, 
see below. 

Relative errors of the elements of our correlation matrices are typically of 0(10"^). This results 
in rather small errors for the fit parameters on a given time interval: masses have relative errors 
of 0{W~^) to 0(10""^), while the components of the energy eigenvectors have errors O(10~^). 
A good fit has to satisfy for our case the following requirements: 

1. The quantum number pattern of the fitting operators has to be as expected, i.e. the 
relative signs of the components of the fitted vectors are the same as those of the rows of 
the maximal twist matrix. 

2. We exclude the results from fit intervals, where the relative errors of the masses are sub- 
stantially higher than the typical errors. With our statistics, this means 1% for the K 
meson mass and 5% for the other masses. Only a few fit intervals turn out to be affected 
by this choice. 

3. The fit ranges [tmin,tmax] are restricted by applying cuts in tmin and {tmax - tmin) such 
that a reasonable "plateau" of the fit values emerges, always keeping a sufficiently large 
number of fit ranges in the sample, typically about 30 to 80. 

After selecting a set of good fits by these criteria a histogram distribution of the fit values has 
been defined by attributing a weight exp(— Xc/dof) to the entries in case of the kaon, and a 
weight l/(Xc/dof) in the other channels. The exponential suppression is in general preferable, 
since it gives robust results but can only be applied for very good fits and plateaus, which is 
the case for the kaon. In order to combine statistical and systematic errors, the entries in the 
distribution were not attributed to a single point but uniformly to the points on the interval 
\Pi — Spi , Pi + Spi] ■ For each final quantity, the quoted value is then the position of the median of 
the resulting distribution. The error is given by a symmetric interval around the median such 
that 68% of the distribution is contained in it. 



We report on single-state, three-state, and four-state fits with a 4 x 4 correlation matrix of 
Gaussian smeared operators. For completeness, we also show the results of three-state fits with 
an 8 X 8 matrix of Gaussian smeared and local operators. All results are summarized in Table [3l 

As shown in table [3] the four-state fit to a 4 x 4 matrix gives one state in each of the channels 
J'^ = 0~ and J = 0^, with both strange and charmed quarks. On the other hand, errors 
are typically larger and/or the light states have higher masses than in the 1x4 and 3x4 fits. 
Therefore, as final results we quote the K meson mass from the 1 x 4 fit and the D meson mass 
from the 3 x 4 fit. 

One can verify a posteriori how well the quantum number content of each fitted vector corre- 
sponds to the expected one. This is simply given by the scalar product of the unit vector in the 
direction of the fitted vector with the row of the matrix in eq. (jlSp that gives the expected vector 
in the continuum limit at maximal twist. For this, we remind that the K meson, strange 0"^ 
state, D meson and charmed 0^ state correspond to the rows 1, 3, 2 and 4, respectively. Table H] 
shows that the fitted vectors are actually well saturated by the expected quantum numbers, 
with scalar products close to 1 in all cases. 
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Ensemble 


M xD 


arriK 


arriK* 


anriD 


amo'^ 


A40.32 


1x4 


0.25542(67) 










3x4 


0.25853(88) 


0.448(13) 


0.903(20) 






4x4 


0.26272(62) 


0.4905(60) 


0.939(46) 


1.09(15) 




3x8 


0.2627(23) 


0.478(18) 


0.885(22) 




B35.32 


1x4 


0.21766(64) 










3x4 


0.21864(51) 


0.422(11) 


0.835(20) 






4x4 


0.2226(69) 


0.449(24) 


0.896(70) 


1.19(10) 




3x8 


0.2203(13) 


0.4369(94) 


0.814(18) 





Table 3: Masses of the K, Kq, D and Dq in lattice units, resulting from the fits to the correlation 
matrices with several eigenstates for the ensembles A40.32 and B35.32. The label M x D means 
a fit with M eigenstates to a -D x D matrix. 



Ensemble 


MxD 


ZK 


Zk* 


ZD 


ZD* 


A40.32 


1x4 


0.98659(6) 










3x4 


0.9871(2) 


0.9896(17) 


0.9392(78) 






4x4 


0.9870(2) 


0.9845(23) 


0.9929(1) 


0.9830(133) 


B35.32 


1x4 


0.98518(8) 










3x4 


0.9847(1) 


0.9772(33) 


0.9518(94) 






4x4 


0.9848(1) 


0.9770(21) 


0.9777(86) 


0.9732(110) 



Table 4: Saturation of the fitted states with the expected quantum numbers, for the four-states 
fits of table El measured by the scalar product z (see text). A value z = 1 indicates complete 
saturation. Mean values and errors for z are determined analogously to masses. 



3.3 Method 3: parity and flavor symmetry restoration 



This third method is a generalization of the "parity restoration method" originally introduced 
for the twisted mass formulation with two degenerate quarks [531 EH EH] . In the N^ = 2 setup the 
twist angle can be determined by requiring that the operators reproducing the correct definition 
of the chiral currents in the continuum limit (physical chiral currents) possess the appropriate 
transformation properties under parity. This condition allows to fix the twist angle for the 
degenerate light quark doublet and the correctly normalized physical currents. We generalize 
the method to the case of bilinear densities with mixed heavy-light flavor composition, used 
here for the determination of the K and D meson masses. A first account of this method can 
be found in [35]. As an outcome, approximations of the physical operators in eq. ([9|) can be 
constructed, from which the masses in the four heavy-flavor and parity channels can be extracted 
by conventional techniques. 



Consider the four-by-four correlation matrix of the renormalized lattice operators in eq. (jM 



C«(t2-tl) 



0^(t2) ^ (O'^ih))^' 



(43) 
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After rewriting the renormalized lattice operators in terms of the bare ones one obtains 
C^{t2-ti) = d[ag( Zp, Zp, Zs, Zs) Cit2 - ti)dmg( Zp, Zp, Zs, Zs) , (44) 



where C(i2 — ^i) is the correlation matrix defined in eq. (j20p . the starting point of the previous 
two methods. The transformation properties of the correlation matrix (|43p can be read from 
eq. p3|) . implying that the correlation matrix of the physical operators Q is given by 

Miu;i,coh)dmg(^Zp,Zp,Zs,Zs)cdmg(^Zp,Zp,Zs,Zs)M^iu;i,coh) , (45) 

where, we recall, the general orthogonal twist rotation matrix Ai{uji,ujh) is given by (J12p . Since 
we are working at maximal twist, we are supposed to insert cji = cjh = it/2 in the rotation matrix 
of eq. (|45p . However, differently from the previous two methods and accounting for the presence 
of 0{a) effects, we treat the two twist angles, along with the renormalization factors Zp and Zs, 
as free parameters. We will return to this point in the following. These free parameters can be 
determined by imposing that the physical operators indeed possess the appropriate parity and 
flavor quantum numbers of their associated channel. This in particular implies that the physical 
correlation matrix of eq. (j45p should be diagonal 



CX)., = 0, j^k. (46) 

Since C{t2—ti) is a symmetric matrix (see section [2.2.ip . the matrix in eq. (|45p is by construction 
symmetric and eq. (I46p actually amounts to only six independent conditions. The latter can be 
rearranged as follows 

(47) 
{+Cn - C22){Zp/Zs) + (-C33 + Cu){Zs/Zp) 



72 

Zp 


C34 


72 ~ 


C12 


ctg{uji) 


\ 
= +- 


ctg(a;h) 


= + 


tan(t<;; + ujh) - 



2(Ci3 — C24:) 

-Cii - C22){Zp/Zs) + (+C33 — Cii){Zs/Zp) 
2(Ci4 — C23) 

Ci4 + C23 + Ci3 + C24 



(+C11 + C22){Zp/Zs) + (-C33 - CiA){Zs/Zp))/2 + Ci2{Zp/Zs) - C^^{ZslZp) 
tan(w/ - ujh) = 



(48) 
(49) 

(50) 



Ci4 + C23 — Ci3 — C24 ,^^. 

^{+Cn + C22){Zp/Zs) + (-C33 - Cu){Zs/Zp))/2 - Ci2{Zp/Zs) + C34(^s/^p) ^ ' 
tan(a;;) ^ C13 + C24 .^2) 

tan{ujh) ~ CM + C23' 

Observe that the right hand sides of (|48p to (j52p are fully determined by the ratio Zp/Zs, i.e. 
they do not depend individually on either Zp or Zs- 
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In Figure [6] we report on the ratios of correlators on the right hand sides of the conditions (j47p 
to (j52p as functions of the time separation t = t2 — ti, for the ensemble A40.32 and the original 
operators without Gaussian smearing. The ratios appear to approach a plateau after a transient: 
from these plateaus we determine the unknown parameters Zp/Zs, uJi and ujh- 

Notice that the time dependence of the ratios is an 0{a) discretization effect, and therefore 
not predicted by eqs. ()47p -( f52|) . which were derived in the continuum limit. For large times the 
lightest eigenstate of the lattice transfer matrix, corresponding to the kaon in the continuum, 
is supposed to saturate the spectral decomposition of the correlation matrix C{t2 — ti) (see 
eqs. (j33|) and ^^)- Assuming a single intermediate state, the six conditions (J17|) - (j52|) are not 
independent any more and in particular the first three of relations ()47p - ()49p are equivalent to 
(|50p - (|52p . Parity and flavor restoration amounts in this case to requiring that the three physical 
operators associated to the heavier channels have no projection on the lightest state, namely 



J](fi|0(f)(x,i)|i^) 



Y,{n\o^;^^\^,t)\K) 



Y,{mo'§^\^,t)\K) 







(53) 



We also observe that this procedure, which relies on asymptotic times, is supposed to be optimal 
from the point of view of the cutoff effects: at large times, contribution from high-mass inter- 
mediate states, which are expected to introduce large discretization effects in the correlator, is 
suppressed. A similar argument was used when tuning the theory to maximal twist in the light 
sector, see [2]. 

We determine Zp/Zs, ^i a-nd uJh by using the relations (l47M9p . while the remaining relations 
serve for cross checking of the results. The latter are reported in Table [5j We observe an 
excellent agreement between the different determinations of the twist angles from (J48M9P and 
(|50ti5ip . respectively, confirming that a single intermediate state contributes. The quality of 
the agreement deteriorates, of course, when the parameters are estimated at smaller temporal 
separations outside the asymptotic region. Notice that the ratio tan{uji) / tan{u}fi) is in all cases 
compatible with zero, since ojh ~ tt/2. Note instead that the value of the light twist angle 



eqs. 


Zp/Zs 


eqs. 


ui/tt 


w/i/vr 


tan(c;j;)/tan(a;/i) 


ensemble A40.32 


m 


0.6575(14) 


iQEEm]) 


0.6504(21) 


0.4980(8) 


-0.012(5) 


m 


same value 


mm 


0.6498(22) 


0.4990(10) 


-0.006(5) 


m 


same value 


(52) 


- 


- 


-0.009(5) 


ensemble B35.32 


(|47D 


0.6793(22) 


iQEEii]) 


0.6453(34) 


0.5005(8) 


0.003(5) 


(|47P 


same value 


dSOlEI]) 


0.6467(29) 


0.5007(9) 


0.005(6) 


m 


same value 


m 


- 


- 


0.005(5) 



Table 5: summary of different determinations of the ratio of renormalization factors and of 
the twist angles with point-like operators (no Gaussian smearing); the first and third column 
indicate the equations used for the determination of the quantities in the corresponding line. 

in Table [5] significantly deviates from the expected value tt/2. In order to understand this 
discrepancy it is useful to recall that the theory is tuned to maximal twist by requiring the 
vanishing of the untwisted PC AC quark mass m^^^ ^ in the light quark sector, see eq. (JH)). 
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Figure 6: ratios of correlators corresponding to the right hand sides of the conditions (j47|) to 
()52p as functions of the temporal separation t for ensemble A40.32 with point-like operators (i.e. 
no Gaussian smearing); the lines indicate the fits in the asymptotic regime. 



This can be shown to be equivalent [55] to requiring parity restoration in the light quark sector. 
One constructs in this case the physical vector current as follows [55] 



ypi^i^) 0^ cos{uji)ZvV^^^+{x) -is{n{i^i)ZAA^^^+{x) , (54) 

where the bilinear of the lattice fields A^'~^{x) is defined in eq. ([5|) and, analogously, 

(55) 



y(0 



X^''h^.X^''^ , 



and Za, Zy are the respective renormalization constants in the massless scheme. The twist 
angle loi is fixed in this case by the condition 



j;(j]|yj')+(x,t)|^) 







(56) 



from which one obtains 



ctg(w;) 



za my I 



PCAC 



X' 







(57) 
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From this we can conclude that our maximal twist condition nv^^ = amounts to ui = '7r/2, if 
the condition (|56|) is assumed. This must be confronted with the conditions (j53|) presently used 
to fix the twist angles ujj^ We conclude that the deviation of uji from 7r/2 should be attributed 
to different 0{a) effects in the pion and kaon sectors. 

We stress that the prescription of eq. (JH), which is based on the charged pion state, is to be 
preferred for tuning the theory to maximal twist, since it ensures the smallest O(a^) discretiza- 
tion errors in physical quantities [32]. Nevertheless, for the determination of the masses in the 
heavy- light meson sector, we use the values of the twist angles obtained from (|47p to (j52p . since 
they deliver optimal projecting operators as defined in eq. (I13p . with the smallest heavy flavor 
and parity violations. The relation in ()57p can also be enforced for the present determination of 
the light twist angle with heavy-light quark bilinears, and the cutoff effects can be absorbed in a 
lattice redeflnition of the PCAC quark mass, m^^/j = 'mF^^'^ + 0{a). For the ensemble A40.32, 

we get for instance Z^in^^) /fJ- ^ —0.5, a pretty large valuqj. The analogous of relation (|57p 
for the heavy twist angle reads 

ctgK) = ^1— . (58) 

The heavy twisted mass /^o- replaces the light twisted mass /i, explaining why lo^ is very close 
to 7r/2: since /^o- ^ fJ-, the non-zero value of rfi^^^'^ only results in a small deviation of uj^ from 
maximal twist. When inserting the above estimate in (j58p we indeed obtain Uh = 0.4956. 

The ratio of normalization factors Zp/Zs and the twist angles loi and ujh allow to determine the 
physical operators up to an overall renormalization (bare physical operators). We choose this 
renormalization to be Zp, so that (cf. eqs. ()13p and ([14 



OtZ'' = Zp'Ol, = M{iOi,ujh)diag(^Zp/Zs,Zp/Zs,hl)o^. (59) 



^ph 



Observe that in the case of the negative parity densities, eq. (|59p corresponds to the conventional 
relation between renormalized and bare operators 

^(h,75)6are ^ Z p' Of^'"'^ ^ ; (60) 

on the other hand, the conventional definition for the bare scalar densities, for which 

^^{h,l)bare,conv. ^—1 ^^(h,l) i? /„^n 

'-^ph - '^S '-^ph \^^) 

holds, is related to the definition (|59p by a finite renormalization 

^^{h,l)bare,conv. y ly ^^{h,l)bare (R')\ 



^In the asymptotic regime, where only the kaon state is considered as intermediate state, the light twist angle 
LJi is fixed by the vanishing of the first two matrix elements in (|53|l : this is so because, in this regime, the two 
conditions can be proven to imply in particular relations H47p and (|48|l (analogously, ujh is fixed in particular by 
the vanishing of the second and third matrix element). 

^For comparison, in the tuning procedure we require ZaI^ti^^) I/m ^ 0.1. 
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Of course, with Zp/Zg at hand both definitions can be computed. 

Figure [7] shows the diagonal and off-diagonal correlators of the bare physical operators for the 
ensemble B35.32 with Gaussian smearing, which are the ones used for the final computation of 
all masses. A general feature is that starting from time separation t > 5 most of the off-diagonal 
elements become small and compatible with zero within statistical errors. An exception is 



the matrix element (Ci^'^^ 



smallest diagonal elements in the (c, 
this observation. 



(^ph )^)' which remains large and comparable in size with the two 



-/+) sectors. At the moment we have no explanation for 



It should also be noted that, following the arguments of |47l I46| . 0{a) improvement can only 
be expected for the diagonal elements of the physical correlation matrix. Since the twist angles 
and the ratio ZpjZg are obtained from conditions on the off-diagonal elements, one should a 
priori expect 0(a) discretization errors for these quantities. However, it should be stressed that 
for physical quantities such as meson masses and decay constants, which are extracted from the 
diagonal matrix elements, 0{a) improvement is at work. The mass of the low- lying state in each 
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Figure 7: Bare physical correlators for the ensemble B35.32 with Gaussian smeared operators. 



of the four different channels can now be extracted by standard techniques from the diagonal 
correlator of the appropriate operator in (I59p . The effective masses for the four channels and 
the two ensembles are reported in Figure [HI for negative parity, and Figure JH for positive parity. 
The final values for all masses are obtained by applying single-mass fits with a cosh function 
in the asymptotic regime. Also in this case the statistical error of the fitting parameters is 
determined by the linearization method of [50]. The starting time tmin for the fits was chosen 
by requiring x^/dof < 1. 

The plateaus for the charmed meson states are generally quite short, since the noise sets in 
early, typically around t > 11. This is, however, expected. For those temporal separations the D 
correlator is only a small fraction of the kaon correlator, as shown in the left panel of Figure [7l 
On the other hand, the D correlator results from a linear combination of the correlators of 
the twisted basis x-fidd bilinears in eq. (jlOp . all dominated by the kaon. This means that the 
condition (|53p can only be fulfilled through a cancellation of large terms, one of the results being 
the comparably small D correlator. The latter inherits the statistical fluctuations of the original 
bilinears and a large relative error is the consequence. As already stated many times in this 
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Figure 8: The effective masses in tlie pseudoscalar channel, with Gaussian smeared operators, 
for the ensembles A40.32 (left) and B35.32 (right). The error bands indicate the total error, 
statistical plus systematic. 



paper, this is an inherent problem in our twisted mass setup, where the D is actually a highly 
excited state in the mixed {s/c, — /+) heavy-light meson sector. 

In the case of the D meson we attempt to estimate the systematic error produced by possible 
residual contributions of excited states and the influence of an unphysical mixing with the rather 
light Kq statqj. We apply a procedure analogous to the one of section 13. 2^ and consider the 
spread of results by including all good fits (those with high significance) obtained by varying 
the fit interval [tminjimax]- The resulting systematic error is much larger than the statistical 
one, and decreases on the finer lattice. This is reflected by the better quality of plateaus for the 
ensemble B35.32, as compared to A40.32, see Figure [8l 

The numerical results for all masses are listed in Table [6l For Kq different plateaus could be 
identified for the effective mass. In this case the value for each plateau is reported. It is unclear 
at this stage, whether this multi plateau behavior reflects the physical structure of QCD states 
in this sector, or is just a statistical effect, as also discussed at the end of section ISTTl 

We conclude the illustration of this method by briefly discussing its generalization to the case of 



*The mixing with the kaon has been ehminated by construction. 
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Figure 9: The effective masses in the scalar channel, with Gaussian smeared operators, for the 
ensembles A40.32 (left) and B35.32 (right). The error bands indicate the total error, statistical 
plus systematic. 



4 X n operators, the immediate application being the one considered in the previous section with 
both local and smeared operators. The obvious route would just be to diagonalize each 4x4 
correlation sub-matrix with homogeneous composition (e.g. local or smeared operators only) as 
we have done so far. As a result, the twist angles and the Zp/Zs factors are obtained for each 
set; observe that the Zp/Zs factors are heavily affected by the smearing, which brings the former 
closer to one. Also the twist angles are expected to differ, due to different 0{a) effects for local 
and smeared operators. Once these parameters are known, the physical correlation matrices 
with mixed local/smeared operators can be reconstructed, too. However, this procedure is not 
expected to be optimal for the latter correlation matrices, since the parameters are adjusted 
to optimize the correlation matrices with homogeneous composition. A better way would be 
to apply an independent diagonalization, with new parameters, of the matrices with mixed 
local/smeared composition. 
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Ensemble 


amx 


^1)^2 


arriK* 


arriD 


aniD* 


A40.32 


0.25668(35) 


7-8 


0.452(8) 


0.909(4)(22) 


1.029(26) 






9-12 


0.431(12) 






14-32 


0.37(5) 


B35.32 


0.21842(33) 


7-10 


0.476(8) 


0.823(4)(14) 


0.968(16) 






11-15 


0.437(23) 






16-32 


0.358(39) 



Table 6: Masses of the K, Kq, D and Dq mesons in lattice units, obtained with the parity and 
flavor restoration method, and using Gaussian smeared operators. The third row contains the 
temporal separations used for the determination of mx* ■ 



4 Conclusions 

We have proposed and compared three methods to determine mx and ttt-d in A^f = 2 -|- 1 -|- 1 
twisted mass lattice QCD. The computation of these masses is less straightforward in this case, 
since parity and flavor are not good quantum numbers. We have therefore explored strategies 
to extract the desired states and have developed three distinct methods all of which exploit 
the exponential fall-off of correlation matrices for suitably chosen heavy-light meson creation 
operators. Method 1 amounts to solving a generalized eigenvalue problem, method 2 is equivalent 
to fitting a linear superposition of exponentials and method 3 transforms the correlators to the 
physical basis by means of the twist rotation. Results for rriK and tti/j obtained with the three 
methods and for both ensembles investigated here are summarized in Table [7] and visualized 
in Figure [TOl Since the kaon is the lightest state in the combined {s/c, —/+) sector, the 





Method 1 


Method 2 


Method 3 


Ensemble A40.32 


arriK 
arriD 


0.2567(2) 
0.922(11) 


0.25554(88) 
0.901(21) 


0.25668(35) 
0.909(22) 


Ensemble B35.32 


arriK 
avfiD 


0.2184(3) 
0.829(8) 


0.21768(84) 
0.835(20) 


0.21842(33) 
0.823(15) 



Table 7: Comparison of the results for rriK and mo obtained with the three methods exposed 
in this work, for both ensembles. 



computation of its mass is rather simple and we obtain precise values for rriK with errors < 0.4% 
including statistical and systematical uncertainties. Moreover, within these errors all three 
methods yield very compatible results which is very reassuring. 

In contrast to ?7i/<, the mass of the D meson is difficult to determine, because in our twisted mass 
setup the D meson is a highly excited state in the combined (s/c, —/+) sector. However, also in 
this case our three methods yield results, which are in excellent agreement within the combined 
statistical and systematical errors, whose relative magnitudes are < 2.5%. Therefore, we are 
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Figure 10: Comparison of the results for rriK (top) and ttid (bottom) obtained with the three 
methods exposed in this work, for both ensembles. The results for methods 1 to 3 are shown 
from left to right. 



confident that we are able to obtain reliable estimates for iriD without resolving all the low lying 
(multi particle) states below the D meson. The latter would require to compute correlation 
matrices with a significantly larger operator basis and with extremely high statistical precision, 
an endeavor, which hardly seems to be feasible. It is therefore very important that already with 
the much smaller correlator matrix employed here, one can obtain a satisfactory estimate of the 
D meson mass. 

The errors we obtain with our three methods differ by factors of around 2 to 4, originating 
from the fact that the three methods estimate the systematic error in different ways. While 
method 2 (fitting exponentials) tends to yield the largest error, its procedure to determine the 
systematic error is also the most conservative: the error is computed from the spread of a 
large set of fit results corresponding to different fitting ranges. In contrast to that method 1 
(solving a generalized eigenvalue problem) estimates the corresponding error by just taking two 
"neighboring fitting ranges" into account. Consequently, the total error is somewhat smaller. 

We stress that as far as K physics is concerned, our analysis shows that this sector can be 
analyzed in the unitary setup without problems. This provides a very good perspective to 
compute corresponding decay constants and also the strange baryon spectrum in the future. For 
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charm physics, the situation is different and it will be quite difficult to extract reliable physics 
results in the charm sector from the unitary setup. Here, we plan to employ a mixed action 
approach by using an Osterwalder-Seiler (OS) valence quark action [M]. This has the advantage 
|56| that there is no flavor mixing and that the valence quarks stay as close as possible to the sea 
twisted mass quarks, e.g. there is no need to re-tune k to realize maximal twist. The idea is to 
match the K and D meson masses between the unitary setup and the valence OS quarks. After 
this matching step further physical quantities such as decay constants will then be computed 
with OS quarks. The matching condition will guarantee that in the continuum limit we recover 
the situation of a unitary setup. Of course, it needs to be seen, whether discretization errors in 
this strategy remain small. Investigations in this direction are in progress. 

With respect to the matching of K and D meson masses between the unitary setup and the 
valence OS quarks, the outcome of our work in this paper is extremely important. The fact that 
we can compute the K meson with high accuracy and the D meson with acceptable precision 
in the unitary setup is a necessary prerequisite to allow for applying such a matching condition. 

Instead of matching the K and D meson masses in the sea and valence sectors, one can directly 
match the renormalized strange and charm quark masses [M]. The latter can be determined in 
the sea sector by using eq. ([6]). Only the finite ratio Zp/Zs is needed as an input for the matching. 
We have shown in this paper one possible way to determine this quantity, which is specific for the 
twisted mass setup. In compliance with the massless quark renormalization scheme, however, 
the extrapolated value of Zp/Zg for four massless quarks is required. We mention here that 
the ETMC has started a dedicated program to evaluate the renormalization constants for our 
A^f = 2 + 1 + 1 setup in the massless quark limit. Once the relevant renormalization constants 
will be available, this information will be used for an alternative tuning of the mass parameters 
in the valence sector. This can result in different values of the valence quark masses with 
respect to the procedure relying on the hadron masses, and hence to different cut-off effects for 
the resulting mixed action theory. Employing both matching conditions can therefore be used 
to have independent computations for physical observables and will provide a most valuable 
cross-check of the way this setup approaches the continuum limit. 



Acknowledgments 

The computer time for this project was made available to us by the John von Neumann-Institute 
for Computing (NIC) on the JUMP, Juropa and Jugene systems in Jiilich and apeNEXT 
system in Zeuthen, BG/P and BG/L in Groningen, by BSC on Mare-Nostrum in Barcelona 
(www.bsc.es), and by the computer resources made available by CNRS on the BlueGene system 
at GENCI-IDRIS Grant 2009-052271 and CCIN2P3 in Lyon. We thank these computer centers 
and their staff for all technical advice and help. 

This work has been supported in part by the DFG Sonderforschungsbereich/ Transregio 
SFB/TR9-03 and the EU Integrated Infrastructure Initiative Hadron Physics (I3HP) under con- 
tract RII3-CT-2004-506078. We also thank the DEISA Consortium (co-funded by the EU, FP6 
project 508830), for support within the DEISA Extreme Computing Initiative (www.deisa.org). 



26 



References 

[1] ETM Collaboration, Ph. Boucaud et al, "Dynamical twisted mass fermions with light 
quarks," Phys. Lett. B 650, 304 (2007) [arXiv:hep-lat/0701012]. 

[2] ETM Collaboration, Ph. Boucaud et al, "Dynamical twisted mass fermions with light 
quarks: simulation and analysis details," Comput. Phys. Commun. 179, 695 (2008) 
[arXiv:0803.0224 [hep-lat]]. 

[3] ETM Collaboration, R. Baron et al, "Light meson physics from maximally twisted mass 
lattice QCD," arXiv:0911.5061 [hep-lat]. 

[4] ETM Collaboration, C. Alexandrou et al, "Light baryon masses with dynamical twisted 
mass fermions," Phys. Rev. D 78, 014509 (2008) [arXiv:0803.3190 [hep-lat]]. 

[5] ETM Collaboration, C. Alexandrou et al, "The low- lying baryon spectrum with two dy- 
namical twisted mass fermions," Phys. Rev. D 80, 114503 (2009) [arXiv:0910.2419 [hep-lat]]. 

[6] ETM Collaboration, B. Blossier et al, "Light quark masses and pseudoscalar decay con- 
stants from A^'f = 2 lattice QCD with twisted mass fermions," JHEP 0804, 020 (2008) 
[arXiv:0709.4574 [hep-lat]]. 

[7] ETM Collaboration, B. Blossier et al, "Pseudoscalar decay constants of kaon and D- 
mesons from A^f = 2 twisted mass lattice QCD," JHEP 0907, 043 (2009) [arXiv:0904.0954 
[hep-lat]]. 

[8] ETM Collaboration, V. Bertone et al, "Kaon oscillations in the Standard Model and 
Beyond using A^f = 2 dynamical quarks," PoS LAT2009 (2009) 258. [arXiv:0910.4838 
[hep-lat]]. 

[9] ETM Collaboration, B. Blossier et al, "A proposal for i?-physics on current lattices," 
arXiv:0909.3187 [hep-lat]. 

[10] ETM Collaboration, B. Blossier et al, "/^ and /b^ with maximally twisted Wilson 
fermions," [arXiv:0911.3757 [hep-lat]]. 

[11] ETM CoUaboration, K. Jansen, C. Michael, A. Shindler and M. Wagner, "The Static-hght 
meson spectrum from twisted mass lattice QCD," JHEP 0812, 058 (2008) [arXiv:0810.1843 
[hep-lat]]. 

[12] ETM Collaboration, C. Michael, A. Shindler and M. Wagner, "The continuum limit of the 
static-light meson spectrum," arXiv: 1004.4235 [hep-lat]. 

[13] ETM Collaboration, B. Blossier, M. Wagner and O. Pene, "Lattice calculation of the 
Isgur-Wise functions ri/2 and T3/2 with dynamical quarks," JHEP 0906, 022 (2009) 
[arXiv:0903.2298 [hep-lat]]. 

[14] ETM Collaboration, R. Frezzotti, V. Lubicz and S. Simula, "Electromagnetic form factor 
of the pion from twisted-mass lattice QCD at A^f = 2," Phys. Rev. D 79, 074506 (2009) 
[arXiv:0812.4042 [hep-lat]]. 



27 



[15] ETM Collaboration, V. Lubicz, F. Mescia, S. Simula et al, "ET — )■ irii' semileptonic form 
factors from two-flavor lattice QCD," Phys. Rev. D80 (2009) 111502. [arXiv:0906.4728 
[hep-lat]]. 

[16] ETM Collaboration, S. Di Vita et al., "Vector and scalar form factors for K- and D-meson 
semileptonic decays from twisted mass fermions with A'^f = 2," PoS LAT2009 (2009) 257. 
[arXiv:0910.4845 [hep-lat]]. 

[17] ETM Collaboration, C. Alexandrou et al, "Nucleon form factors with dynamical twisted 
mass fermions," PoS LATTICE2008, 139 (2008) [arXiv:0811.0724 [hep-lat]]. 

[18] ETM Collaboration, R. Baron, S. Capitani, J. Carbonell, K. Jansen, Z. Liu, O. Pene 
and C. Urbach, "Moments of meson distribution functions with dynamical twisted mass 
fermions," PoS LAT2007, 153 (2007) [arXiv:0710.1580 [hep-lat]]. 

[19] ETM Collaboration, C. Michael and C. Urbach, "Neutral mesons and disconnected dia- 
grams in twisted mass QCD," PoS LAT2007, 122 (2007) [arXiv:0709.4564 [hep-lat]]. 

[20] ETM Collaboration, K. Jansen, C. Michael and C. Urbach, "The rj' meson from lattice 
QCD," Eur. Phys. J. C 58, 261 (2008) [arXiv:0804.3871 [hep-lat]]. 

[21] ETM Collaboration, C. McNeile, C. Michael and C. Urbach, "The to-p meson mass splitting 
and mixing from lattice QCD," Phys. Lett. B 674, 286 (2009) [arXiv:0902.3897 [hep-lat]]. 

[22] ETM Collaboration, D. B. Renner and X. Feng, "Hadronic contribution to ^r — 2 from 
twisted mass fermions," PoS LATTICE2008, 129 (2008) [arXiv: 0902. 2796 [hep-lat]]. 

[23] ETM Collaboration, X. Feng, K. Jansen and D. B. Renner, "The vr"*" tt~^ scattering length 
from maximally twisted mass lattice QCD," Phys. Lett. B 684, 268 (2010) [arXiv:0909.3255 
[hep-lat]]. 

[24] ETM Collaboration, X. Feng, K. Jansen and D. B. Renner, "Scattering from finite size 
methods in lattice QCD," arXiv:0910.4871 [hep-lat]. 

[25] ETM Collaboration, M. Constantinou et al, "Non-perturbative renormalization of quark 
bilinear operators with N{ = 2 (tmQCD) Wilson fermions and the tree-level improved gauge 
action," arXiv:1004.1115 [hep-lat]. 

[26] K. Cichy, J. Gonzalez Lopez, K. Jansen, A. Kujawa and A. Shindler, "Twisted mass, overlap 
and Creutz fermions: cut-off effects at tree-level of perturbation theory," Nucl. Phys. B 800, 
94 (2008) [arXiv:0802.3637 [hep-lat]]. 

[27] XjF Collaboration, K. Jansen, M. Papinutto, A. Shindler, C. Urbach and I. Wetzorke, 
"Light quarks with twisted mass fermions," Phys. Lett. B 619, 184 (2005) [arXiv:hep- 
lat/0503031]. 

[28] ^L^ Collaboration, K. Jansen, M. Papinutto, A. Shindler, C. Urbach and I. Wetzorke, 
"Quenched scaling of Wilson twisted mass fermions," JHEP 0509, 071 (2005) [arXiv:hep- 
lat/0507010]. 

A. M. Abdel-Rehim, R. Lewis and R. M. Woloshyn, "Spectrum of quenched twisted mass 
lattice QCD at maximal twist," Phys. Rev. D 71, 094505 (2005) [arXiv:hep-lat/0503007] . 

28 



[30] ETM Collaboration, C. Urbacli, "Lattice QCD with two light Wilson quarks and maximally 
twisted mass," PoS LAT2007 (2007) 022 [arXiv:0710.1517 [hep-lat]]. 

[31] ETM Collaboration, P. Dimopoulos, R. Frezzotti, G. Herdoiza, C. Urbach and U. Wenger, 
"Scaling and low energy constants in lattice QCD with Nf = 2 maximally twisted Wilson 
quarks," PoS LAT2007 (2007) 102 [arXiv:0710.2498 [hep-lat]]. 

[32] R. Frezzotti, G. Martinelli, M. Papinutto and G. C. Rossi, "Reducing cutoff effects in 
maximally twisted lattice QCD close to the chiral limit," JHEP 0604 038 (2006) [arXiv:hep- 
lat/0503034]. 

[33] ETM Collaboration, P. Dimopoulos, R. Frezzotti, C. Michael, G. C. Rossi and C. Urbach, 
"©(a^) cutoff effects in lattice Wilson fermion simulations," Phys. Rev. D 81, 034509 (2010) 
[arXiv:0908.0451 [hep-lat]]. 

[34] R. Frezzotti and G. C. Rossi, "Chirally improving Wilson fermions. II: Four-quark opera- 
tors," JHEP 0410 (2004) 070 [arXiv:hep-lat/0407002]. 

[35] T. Chiarappa et al, "Numerical simulation of QCD with u, d, s and c quarks in the twisted- 
mass Wilson formulation," Eur. Phys. J. C 50, 373 (2007) [arXiv:hep-lat/0606011]. 

[36] ETM Cohaboration, R. Baron et al, "Status of ETMC simulations with Nf = 2 + 1 + 1 
twisted mass fermions," PoS LATTICE2008, 094 (2008) [arXiv:0810.3807 [hep-lat]]. 

[37] ETM Collaboration, R. Baron et al, "First results of ETMC simulations with Nf = 2-M-M 
maximally twisted mass fermions," arXiv:0911.5244 [hep-lat]. 

[38] ETM Collaboration, R. Baron et a/., "Light hadrons from lattice QCD with light {u,d), 
strange and charm dynamical quarks", arXiv: 1004. 5284 [hep-lat]. 

[39] MILC Collaboration, A. Bazavov et a/., "HISQ action in dynamical simulations," PoS 
LATTICE2008 (2008) 033 [arXiv:0903.0874 [hep-lat]]. 

[40] MILC Collaboration, A. Bazavov et al, "Progress on four flavor QCD with the HISQ 
action," PoS LAT2009 (2009) 123 [arXiv:0911.0869 [hep-lat]]. 

[41] MILC Collaboration, A. Bazavov et al, "Scaling studies of QCD with the dynamical HISQ 
action," [arXiv: 1004.0342 [hep-lat]]. 

[42] C. Pena, S. Sint and A. Vladikas, "Twisted mass QCD and lattice approaches to the Delta(I) 
= 1/2 rule," JHEP 0409 (2004) 069 [arXiv:hep-lat/0405028]. 

[43] A. M. Abdel-Rehim, R. Lewis, R. M. Woloshyn and J. M. S. Wu, "Strange quarks 
in quenched twisted mass lattice QCD," Phys. Rev. D 74 (2006) 014507 [arXiv:hep- 
lat/0601036]. 

[44] Y. Iwasaki, "Renormalization group analysis of lattice theories and improved lattice action: 
two-dimensional non-linear 0{N) sigma model," Nucl. Phys. B 258, 141 (1985). 

[45] ALPHA Cohaboration, R. Frezzotti, P. A. Grassi, S. Sint and P. Weisz, "Lattice QCD 
with a chirally twisted mass term," JHEP 0108, 058 (2001) [arXiv:hep-lat/0101001]. 



29 



[46] R. Frezzotti and G. C. Rossi, "Twisted-mass lattice QCD with mass non-degenerate 
quarks," Nucl. Phys. Proc. Suppl. 128 (2004) 193 [arXiv:hep-lat/0311008]. 

[47] R. Frezzotti and G. C. Rossi, "Chirally improving Wilson fermions. I: 0{a) improvement," 
JHEP 0408, 007 (2004) [arXiv:hep-lat/0306014]. 

[48] C. Amsler et al. [Particle Data Group], Phys. LettB 667, 1 (2008) and 2009 partial update 
for the 2010 edition. 

[49] ALPHA Cohaboration, B. Blossier, M. Delia Morte, G. von Hippel, T. Mendes and 
R. Sommer, "On the generalized eigenvalue method for energies and matrix elements in 
lattice field theory," JHEP 0904, 094 (2009) [arXiv:0902.1265 [hep-lat]]. 

[50] ALPHA Collaboration, U. Wolff, "Monte Carlo errors with less errors," Comput. Phys. 
Commun. 156, 143 (2004) [Erratum-ibid. 176, 383 (2007)] [arXiv:hep-lat/0306017]. 

[51] C. Michael, "Fitting correlated data," Phys. Rev. D 49, 2616 (1994) [arXiv:hep- 
lat/9310026]. 

[52] C. Michael and A. McKerrell, "Fitting Correlated Hadron Mass Spectrum Data," Phys. 
Rev. D 51, 3745 (1995) [arXiv:hep-lat/9412087]. 

[53] F. Farchioni et al., "Exploring the phase structure of lattice QCD with twisted mass 
quarks," Nucl. Phys. Proc. Suppl. 140, 240 (2005) [arXiv:hep-lat/0409098]. 

[54] F. Farchioni et al., "The phase structure of lattice QCD with Wilson quarks and renormal- 
ization group improved gluons," Eur. Phys. J. C 42, 73 (2005) [arXiv:hep-lat/0410031]. 

[55] F. Farchioni et al., "Numerical simulations with two flavours of twisted-mass Wilson quarks 
and DBW2 gauge action," Eur. Phys. J. C 47, 453 (2006) [arXiv:hep-lat/0512017]. 

[56] R. Frezzotti, private communication (2008). 



30 



